{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further Hypothesis Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──\n",
      "✔ ggplot2 3.3.0     ✔ purrr   0.3.4\n",
      "✔ tibble  3.0.1     ✔ dplyr   0.8.5\n",
      "✔ tidyr   1.1.0     ✔ stringr 1.4.0\n",
      "✔ readr   1.3.1     ✔ forcats 0.4.0\n",
      "── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──\n",
      "✖ dplyr::filter() masks stats::filter()\n",
      "✖ dplyr::lag()    masks stats::lag()\n"
     ]
    }
   ],
   "source": [
    "# Select this cell and type Ctrl-Enter to execute the code below.\n",
    "\n",
    "library(tidyverse)\n",
    "\n",
    "set_plot_dimensions <- function(width_choice, height_choice) {\n",
    "    options(repr.plot.width=width_choice, repr.plot.height=height_choice)\n",
    "}\n",
    "\n",
    "cbPal <- c(\"#E69F00\", \"#56B4E9\", \"#009E73\", \"#F0E442\", \"#CC79A7\", \"#0072B2\", \"#D55E00\")\n",
    "\n",
    "set_plot_dimensions(5, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You should see \"Attaching packages\" and some ticks by the packages loaded.\n",
    "# The \"Conflicts\" aren't a problem.\n",
    "\n",
    "# Other problems loading the library? Try running this cell.\n",
    "\n",
    "install.packages('tidyverse')\n",
    "\n",
    "library(tidyverse)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 - Comparing distributions over a categorical variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Parsed with column specification:\n",
      "cols(\n",
      "  temperature = col_double(),\n",
      "  luminosity = col_double(),\n",
      "  radius = col_double(),\n",
      "  spectral_class = col_character(),\n",
      "  type = col_double()\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "# Run this cell to load the data.\n",
    "\n",
    "data <- read_csv(\"stars.csv\")\n",
    "\n",
    "type_key <- c('Brown Dwarf', 'Red Dwarf', 'White Dwarf', 'Main Sequence', 'Supergiant','Hypergiant')\n",
    "spectral_classes <- c('O','B','A','F','G','K','M')\n",
    "\n",
    "data$type <- factor(data$type)\n",
    "data$spectral_class <- factor(data$spectral_class, levels=spectral_classes)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your undergraduate project student, Tunde, is looking at the distributions of spectral class (i.e. colour) for various types of star. \n",
    "\n",
    "He says that the bar charts for white dwarves (type 2) and main sequence stars (type 3) seem similar, but you are not so sure. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAHgCAYAAAB5FxNZAAAEGWlDQ1BrQ0dDb2xvclNwYWNl\nR2VuZXJpY1JHQgAAOI2NVV1oHFUUPrtzZyMkzlNsNIV0qD8NJQ2TVjShtLp/3d02bpZJNtoi\n6GT27s6Yyc44M7v9oU9FUHwx6psUxL+3gCAo9Q/bPrQvlQol2tQgKD60+INQ6Ium65k7M5lp\nurHeZe58853vnnvuuWfvBei5qliWkRQBFpquLRcy4nOHj4g9K5CEh6AXBqFXUR0rXalMAjZP\nC3e1W99Dwntf2dXd/p+tt0YdFSBxH2Kz5qgLiI8B8KdVy3YBevqRHz/qWh72Yui3MUDEL3q4\n4WPXw3M+fo1pZuQs4tOIBVVTaoiXEI/MxfhGDPsxsNZfoE1q66ro5aJim3XdoLFw72H+n23B\naIXzbcOnz5mfPoTvYVz7KzUl5+FRxEuqkp9G/Ajia219thzg25abkRE/BpDc3pqvphHvRFys\n2weqvp+krbWKIX7nhDbzLOItiM8358pTwdirqpPFnMF2xLc1WvLyOwTAibpbmvHHcvttU57y\n5+XqNZrLe3lE/Pq8eUj2fXKfOe3pfOjzhJYtB/yll5SDFcSDiH+hRkH25+L+sdxKEAMZahrl\nSX8ukqMOWy/jXW2m6M9LDBc31B9LFuv6gVKg/0Szi3KAr1kGq1GMjU/aLbnq6/lRxc4XfJ98\nhTargX++DbMJBSiYMIe9Ck1YAxFkKEAG3xbYaKmDDgYyFK0UGYpfoWYXG+fAPPI6tJnNwb7C\nlP7IyF+D+bjOtCpkhz6CFrIa/I6sFtNl8auFXGMTP34sNwI/JhkgEtmDz14ySfaRcTIBInmK\nPE32kxyyE2Tv+thKbEVePDfW/byMM1Kmm0XdObS7oGD/MypMXFPXrCwOtoYjyyn7BV29/MZf\nsVzpLDdRtuIZnbpXzvlf+ev8MvYr/Gqk4H/kV/G3csdazLuyTMPsbFhzd1UabQbjFvDRmcWJ\nxR3zcfHkVw9GfpbJmeev9F08WW8uDkaslwX6avlWGU6NRKz0g/SHtCy9J30o/ca9zX3Kfc19\nzn3BXQKRO8ud477hLnAfc1/G9mrzGlrfexZ5GLdn6ZZrrEohI2wVHhZywjbhUWEy8icMCGNC\nUdiBlq3r+xafL549HQ5jH+an+1y+LlYBifuxAvRN/lVVVOlwlCkdVm9NOL5BE4wkQ2SMlDZU\n97hX86EilU/lUmkQUztTE6mx1EEPh7OmdqBtAvv8HdWpbrJS6tJj3n0CWdM6busNzRV3S9KT\nYhqvNiqWmuroiKgYhshMjmhTh9ptWhsF7970j/SbMrsPE1suR5z7DMC+P/Hs+y7ijrQAlhyA\ngccjbhjPygfeBTjzhNqy28EdkUh8C+DU9+z2v/oyeH791OncxHOs5y2AtTc7nb/f73TWPkD/\nqwBnjX8BoJ98VQNcC+8AAEAASURBVHgB7d0JfBRF2vjxJ3dISEK4NICAeCCC9yqioCKCeIGu\nuru8qOC1oqCIguK5Lx54oOIFHqCsLqIIrseq6K54oQvei3KKyiU3yBUIEJL5z1P79vxnkoEw\nJNPd6frV55Nkuqenq+pbk55nqqurU0LhJCQEEEAAAQQQQACBGhNIrbE9sSMEEEAAAQQQQAAB\nI0CAxRsBAQQQQAABBBCoYQECrBoGZXcIIIAAAggggAABFu8BBBBAAAEEEECghgUIsGoYlN0h\ngAACCCCAAAIEWLwHEEAAAQQQQACBGhZIr+H97dXuli9fLtOmTZO0tDQ54YQTpEmTJpH9bN68\nWaZPnx5Zdh507txZMjIynEX+IoAAAggggAACvhFI8XoerDvuuEO++OIL6dSpkyxcuFAWL14s\n99xzj3To0MEgffbZZ3L77bdLw4YNY9DGjRsneXl5MetYQAABBBBAAAEE/CDgaQ/W/Pnz5dNP\nP5VJkyZJ48aNjcewYcPk8ccfjwRYCxYskLZt28qoUaMS9tq0aZNs27Yt4dcl8wV169aV9PR0\n2bBhQzKz8eW+8/PzTbm0XWxL9erVk507d0pxcbFVVU9NTRVt9x07dsjWrVutqrv+n+v/ux6D\n/HYcSnZDZGZmSk5OjmzZskVKS0uTnZ2v9p+dnS36o2dfysrKfFU253PWV4UKcGE8DbDWr18v\nl19+eSS4UuejjjpKPv74Y9EJ5lNSUkQDrNatW+9VE5SXl/vuDa4HXT34+O0fb6+AE3yRntLV\ndrWx7trmNtZd/4e17hpc2tbuzv/69u3brau70+4aVNvW7vqlwtb/9wQ/EgK/uacB1vHHHy/6\nE52mTp0qbdq0McGVrtcAKysrS4YOHSrz5s0zzw0YMECaNm0a/TJZtWqVjB49OmZd165dTcAW\ns9LjBR1npsnpzfG4OK5mrwdd/bGx7gqtH7i21V3bW5MG17bV3flf1+OXfujalPS9rkl7cmwb\nK+vUNzc3V/RLPsleAU8DrIrsEydOlJkzZ8ozzzxjntIu1pUrV8q+++4rvXr1ko4dO8rkyZOl\nf//+Mn78eNP97uxDT7m98sorzqL5e/DBB5vXxKz0yYL+89mabK27fuDaWnf90HE+eGx732tv\nhv7YmDTAsjXVqVPH1qpT7/8T8HyQu9MSzz//vLz00kty1113yYknnmhW6ykV7ZmqX79+5AD1\nyy+/SJ8+fWTw4MHSs2dP5+Wi3fD6XHTSDzO/vcl1LI5+0KxZsya6qFY8btCggTlN9ttvv1lR\n3+hKNmrUyIxD2rhxY/TqwD/WoFL/f0tKSqwbf6ZBVUFBgRmHZNv4Mw2s9CIkHW+px2abkn7u\n6PgzHQKjp8b9lPQ4RHJPwPMeLO1Cffjhh+WDDz6Qhx56KOaUnp5e0N6r6NSqVSvRN8mKFSui\nV5vTiHpqMTppr5Ye2P2UNGjU5Ld/PDeMtO76Y2Pd1Ze6++vDJtnveecUoR7jbHvPO6fGbK67\njj2zrd2T/T9V2/bv+cCAu+++28xz9dRTT8UEVwq5aNEi01u1dOnSiKsGVtr7U3EMVmQDHiCA\nAAIIIIAAAh4LeBpgTZkyxfRc9e3b11zSquOvnB+N/lu2bGkGST799NOmu1WDKx3IXlhYKF26\ndPGYjuwRQAABBBBAAIH4Ap6eItQB65pGjBhRqXTvv/++OY89aNAgMy7rvPPOM9voKcInn3zS\nPFfpRaxAAAEEEEAAAQR8IOBpgPXcc89VSXDIIYfIhAkTZO3atWZwuA4aJSGAAAIIIIAAAn4W\n8DTASgSm4q1yEnkt2yKAAAIIIIAAAm4KeDoGy82KkhcCCCCAAAIIIOCWAAGWW9LkgwACCCCA\nAALWCBBgWdPUVBQBBBBAAAEE3BIgwHJLmnwQQAABBBBAwBoBAixrmpqKIoAAAggggIBbAgRY\nbkmTDwIIIIAAAghYI0CAZU1TU1EEEEAAAQQQcEuAAMstafJBAAEEEEAAAWsEas1Eo9a0iMsV\nPWzira7muHLAGFfzIzMEEEAAAQS8EKAHywt18kQAAQQQQACBQAsQYAW6eakcAggggAACCHgh\nQIDlhTp5IoAAAggggECgBQiwAt28VA4BBBBAAAEEvBAgwPJCnTwRQAABBBBAINACBFiBbl4q\nhwACCCCAAAJeCBBgeaFOnggggAACCCAQaAECrEA3L5VDAAEEEEAAAS8ECLC8UCdPBBBAAAEE\nEAi0AAFWoJuXyiGAAAIIIICAFwIEWF6okycCCCCAAAIIBFqAACvQzUvlEEAAAQQQQMALAQIs\nL9TJEwEEEEAAAQQCLUCAFejmpXIIIIAAAggg4IUAAZYX6uSJAAIIIIAAAoEWIMAKdPNSOQQQ\nQAABBBDwQoAAywt18kQAAQQQQACBQAsQYAW6eakcAggggAACCHghQIDlhTp5IoAAAggggECg\nBQiwAt28VA4BBBBAAAEEvBAgwPJCnTwRQAABBBBAINACBFiBbl4qhwACCCCAAAJeCKR7kalb\neaalpUl+fr5b2e1RPlomTX4r1x4VvgY2Sk1Ntbbu6enp1tU9JSXFvGsyMjKsq7vzv56VlSX6\nvrcp6XtdU3Z2tmjb25Sc+ubm5kp5eblNVaeuFQQCHWDpm3vnzp0Vquztoh5sNe3YscPbgniU\neygUsrLuzsHWtnZ3Agv9X7St7hpkaIBRVlZmXd2dwNrGujuBdWlpqWl7jw61cbPV4xDJPYFA\nB1j6Yb5t2zb3NPcgJ+cN7rdy7UHRa2QTP7ZJjVRsD3aiQYZt7e70ZOgHrW11d75M6Zc82+ru\nBNYaZNhWd6cHS79QaP1J9grY1W9tbztTcwQQQAABBBBwUYAAy0VsskIAAQQQQAABOwQIsOxo\nZ2qJAAIIIIAAAi4KEGC5iE1WCCCAAAIIIGCHAAGWHe1MLRFAAAEEEEDARQECLBexyQoBBBBA\nAAEE7BAgwLKjnaklAggggAACCLgoQIDlIjZZIYAAAggggIAdAgRYdrQztUQAAQQQQAABFwUI\nsFzEJisEEEAAAQQQsEOAAMuOdqaWCCCAAAIIIOCiAAGWi9hkhQACCCCAAAJ2CBBg2dHO1BIB\nBBBAAAEEXBQgwHIRm6wQQAABBBBAwA4BAiw72plaIoAAAggggICLAgRYLmKTFQIIIIAAAgjY\nIUCAZUc7U0sEEEAAAQQQcFGAAMtFbLJCAAEEEEAAATsECLDsaGdqiQACCCCAAAIuChBguYhN\nVggggAACCCBghwABlh3tTC0RQAABBBBAwEUBAiwXsckKAQQQQAABBOwQIMCyo52pJQIIIIAA\nAgi4KECA5SI2WSGAAAIIIICAHQIEWHa0M7VEAAEEEEAAARcFCLBcxCYrBBBAAAEEELBDgADL\njnamlggggAACCCDgogABlovYZIUAAggggAACdggQYNnRztQSAQQQQAABBFwUIMByEZusEEAA\nAQQQQMAOAQIsO9qZWiKAAAIIIICAiwIEWC5ikxUCCCCAAAII2CGQ7odqLl++XKZNmyZpaWly\nwgknSJMmTWKKtXnzZvn8889F/7Zv316aN28e8zwLCCCAAAIIIICAnwQ878G64447pG/fvvLj\njz/Ku+++KxdffLFMnz49YrRw4ULp2bOnTJ48WWbNmiWXXXaZzJgxI/I8DxBAAAEEEEAAAb8J\neNqDNX/+fPn0009l0qRJ0rhxY2MzbNgwefzxx6VDhw5m+b777pMePXrIwIEDJSUlRV544QUZ\nOXKkvPLKK2bZb6CUBwEEEEAAAQQQ8LQHa/369XL55ZdHgittjqOOOkpWrlwpoVBI1q1bJ3Pn\nzjU9WBpcaTr77LNFTynOmTPHLPMLAQQQQAABBBDwm4CnPVjHH3+86E90mjp1qrRp08b0Tmmg\npSl6TFaDBg0kMzNTVq9eLW3bto28VIO1d955J7KsDw477DBp2bJlzDqvF1JT/xvT5uTkeF0U\nT/LX+ttadx1jaFvdnfd7enq6dXXXOmvKyMiwru56jNakf50vx2aFBb+0vTVlZ2ebtregylRx\nFwKeBlgVyzRx4kSZOXOmPPPMM+apFStWSFZWlvmJ3jYvL080oIpOGnDdfffd0avkzjvvlCOO\nOCJmnV8WCgoK/FIUV8uhH7i21l0/cG2tu37QOh+6rr7hfJCZftDqj43Jti8U0W1ct27d6EUe\nWyjgmwDr+eefl5deeknuvfdead26tWkK/Sawc+fOSs1SVlZW6Ruh9nI99NBDMdsecMABlQKx\nmA08WNDgUD9oKwaIHhTFkyy17TZt2uRJ3l5mWlhYKKWlpVJcXOxlMVzP2wmot2/fLlu3bnU9\nfy8z1P9z/X8vKSmRbdu2eVkU1/PWL8YaXG3ZskV27Njhev5eZlinTh0TUOtxTo93fkp6HCK5\nJ+B5gFVeXi4PP/ywfPDBByZA0jFYTmrYsKF5g+qBOfqbkL5xi4qKnM3MXz2QnXPOOTHrNmzY\nYA5uMSs9XsjNzTUlsO2A67Dr2Dpb667vddvq7pwm0w8a2+quQYYm/ZJoW92dU8P6pcK2ujun\nCDWw1PqT7BXwdJC7sutpPZ2W4amnnjID3KObolmzZqa3Z/bs2ZHVOuhdP6iix2VFnuQBAggg\ngAACCCDgAwFPe7CmTJlieq6GDBliJhHV8VdOateunRmv0q1bNxk3bpwZ+K7fhseOHSvdu3eX\nRo0aOZvyFwEEEEAAAQQQ8JWApwGWTh6qacSIEZVQ3n//fXNasF+/fqJzY+npP+1y10Hr1157\nbaXtWYEAAggggAACCPhFwNMA67nnnqvSQQflPfroo2ZgtF7m7oxhqvKFbIAAAggggAACCHgk\n4GmAlUid8/PzE9mcbRFAAAEEEEAAAc8EPB/k7lnNyRgBBBBAAAEEEEiSAAFWkmDZLQIIIIAA\nAgjYK0CAZW/bU3MEEEAAAQQQSJIAAVaSYNktAggggAACCNgrQIBlb9tTcwQQQAABBBBIkgAB\nVpJg2S0CCCCAAAII2CtAgGVv21NzBBBAAAEEEEiSAAFWkmDZLQIIIIAAAgjYK0CAZW/bU3ME\nEEAAAQQQSJIAAVaSYNktAggggAACCNgrQIBlb9tTcwQQQAABBBBIkgABVpJg2S0CCCCAAAII\n2CtAgGVv21NzBBBAAAEEEEiSAAFWkmDZLQIIIIAAAgjYK0CAZW/bU3MEEEAAAQQQSJIAAVaS\nYNktAggggAACCNgrQIBlb9tTcwQQQAABBBBIkgABVpJg2S0CCCCAAAII2CtAgGVv21NzBBBA\nAAEEEEiSAAFWkmDZLQIIIIAAAgjYK0CAZW/bU3MEEEAAAQQQSJIAAVaSYNktAggggAACCNgr\nQIBlb9tTcwQQQAABBBBIkgABVpJg2S0CCCCAAAII2CtAgGVv21NzBBBAAAEEEEiSAAFWkmDZ\nLQIIIIAAAgjYK0CAZW/bU3MEEEAAAQQQSJIAAVaSYNktAggggAACCNgrQIBlb9tTcwQQQAAB\nBBBIkgABVpJg2S0CCCCAAAII2CtAgGVv21NzBBBAAAEEEEiSQHqS9uuL3aampkqdOnV8URan\nEFomTX4rl1O+ZP9NSUmxtu5+fD8mu72d93taWpp17Z6e/t/Dq/617f89IyPDvLWcv8l+n/lp\n/067Z2VlifPYT+WjLO4JBDrAUkY9sPspaYChyW/lcstI60/d3dL2Ph8nwNK/trU7dRexud1t\ne797f7TxXwkCHWCVl5dLSUmJr9T1W43+4xUXF/uqXG4VRtvExrrn5eVJWVmZdXXXb/C5ublS\nWlpqXd31fz0nJ0d27NhhXd213tnZ2bJ9+3bfHYOTfazTL5GZmZmydetW875Pdn6J7F+PQyT3\nBBiD5Z41OSGAAAIIIICAJQIEWJY0NNVEAAEEEEAAAfcECLDcsyYnBBBAAAEEELBEgADLkoam\nmggggAACCCDgngABlnvW5IQAAggggAAClggQYFnS0FQTAQQQQAABBNwTIMByz5qcEEAAAQQQ\nQMASAQIsSxqaaiKAAAIIIICAewIEWO5ZkxMCCCCAAAIIWCJAgGVJQ1NNBBBAAAEEEHBPgADL\nPWtyQgABBBBAAAFLBAiwLGloqokAAggggAAC7gkQYLlnTU4IIIAAAgggYIkAAZYlDU01EUAA\nAQQQQMA9AQIs96zJCQEEEEAAAQQsESDAsqShqSYCCCCAAAIIuCdAgOWeNTkhgAACCCCAgCUC\nBFiWNDTVRAABBBBAAAH3BAiw3LMmJwQQQAABBBCwRIAAy5KGppoIIIAAAggg4J4AAZZ71uSE\nAAIIIIAAApYIEGBZ0tBUEwEEEEAAAQTcEyDAcs+anBBAAAEEEEDAEgECLEsammoigAACCCCA\ngHsCBFjuWZMTAggggAACCFgiQIBlSUNTTQQQQAABBBBwT4AAyz1rckIAAQQQQAABSwQIsCxp\naKqJAAIIIIAAAu4JEGC5Z01OCCCAAAIIIGCJAAGWJQ1NNRFAAAEEEEDAPQECLPesyQkBBBBA\nAAEELBEgwLKkoakmAggggAACCLgnkO5eVrvPqaysTMaPHy/nnXee5OfnRzbevHmzTJ8+PbLs\nPOjcubNkZGQ4i/xFAAEEEEAAAQR8I+CbAGv06NHy6quvSteuXWMCrJkzZ8rw4cOlYcOGMWgd\nOnQgwIoRYQEBBBBAAAEE/CLgeYC1atUqeeihh+Tbb7+Na7JgwQJp27atjBo1Ku7zrEQAAQQQ\nQAABBPwm4PkYrPvvv19CoZA88MADcW00wGrdunXc51iJAAIIIIAAAgj4UcDzHqyhQ4fKPvvs\nI4sXL47rowFWVlaW6Hbz5s2TNm3ayIABA6Rp06Yx22/dulX0dGJ0aty4sdSrVy96leePU1JS\nTBkyMzM9L4sXBdD621r31NRU6+qelpZm3mY21j09/b+HVzWw7T3vtLvNddcxws7x3otjLXl6\nL+B5gKXB1a6SDnBfuXKl7LvvvtKrVy/p2LGjTJ48Wfr3728GxNetWzfy0qVLl0rfvn0jy/rg\nzjvvlN69e8es88tCgwYN/FIUV8uhB1xb664HXFvrnp2dLfpjY8rJyRH9sTHl5eXZWG1T54KC\nAmvrTsX/K+B5gLW7htAAatKkSVK/fv3IN8BDDz1U+vTpI1OnTpWePXtGXq4fXP369Yss64OD\nDjpINEjzU9IDrQYZfiuXW0bl5eWyZcsWt7LzTT76QbNz504pKSnxTZncKIj2XOXm5kppaals\n27bNjSx9k4f+n+v/+/bt22XHjh2+KZcbBdEvExpQ6/td3/c2JT3joj2WelZFr473U7I54PWi\nHXwdYGn3qvZeRadWrVpJo0aNZMWKFdGrzVWGgwYNilm3YcMGKS4ujlnn9YL+8+mB12/lcstF\nAywb664HNj3Y2lZ3PU3mBFi21V3/1zXA0uDKtrprvTXA0uDSti8VzjAIDbD0i4WfEgGWu63h\n+SD33VV30aJFprdKT/85SQOrNWvWVBqD5TzPXwQQQAABBBBAwGsBXwdYLVu2NN+Cnn76aVm/\nfr3ptdL5sgoLC6VLly5e25E/AggggAACCCAQV8DXAZaWWE/7LVy40MzwrgPdly1bJk8++aS1\ng0bjtiIrEUAAAQQQQMBXAr4Zg9WiRQuZNm1aJZxDDjlEJkyYIGvXrjUzt3NlRiUiViCAAAII\nIICAzwR8E2BV5VLxVjlVbc/zCCCAAAIIIICAVwK+P0XoFQz5IoAAAggggAACeytAgLW3crwO\nAQQQQAABBHYpoNNUPPjgg+bK/11uFOAnCLAC3LhUDQEEEEAAAa8ERowYITfffLN1E+063gRY\njgR/EUAAAQQQQKDGBGybxb8iXK0Z5F6x4CwjgAACCCCAgD8FJk6cKB9++KEp3EMPPSRHHXWU\nnHvuufLwww9HHkeXXCcR1zkvzz77bKlXr5653/A111xj9vHPf/7TzH951llnySmnnGLuhhL9\n2h9++EFeffVVmTt3rjRv3tzs49RTT43exJPH9GB5wk6mCCCAAAIIBFdgzpw5snjxYlPBL7/8\nUmbNmiX5+fny+uuvy3XXXSehUCim8n/961/l7rvvFr2v8IIFC+Suu+6SSy+9VK699lqznb7+\ntNNOkyFDhsS87plnnpHf/e53ZqyX3pbqk08+MRORV9wu5kUuLRBguQRNNggggAACCNgiMGzY\nMBMgaX21d0kHu2vSoElvf6eBUHR68cUXpVOnTqL3G3bS119/LV999ZWMGzdOtBdr+PDhMnLk\nSBOk6TY//fSTCdZOPPFEMwn5W2+9Jd98843cdtttor1mU6dOdXblyV8CLE/YyRQBBBBAAAH7\nBHr37m0mDR8/fnyk8hpEzZs3T/r27RtZpw+uv/562X///SPrbrzxRtE5MV977TWzTm+dp71W\nQ4cONeudDXW7jIwMGTVqlLPKk7+MwfKEnUwRQAABBBCwT6Bx48aiY6kmT55sbnuXnZ0t2nuV\nm5srF154YQxI+/btY5YzMzNF7+7y3XffmfXz58+XlJQUefbZZ2Xs2LEx2+bk5MiPP/4Ys87t\nBXqw3BYnPwQQQAABBCwW0NOEGzdulLffflt0rqxXXnlFzj//fKlbt26Mig52r5g0cNLXatJb\n6GVlZUl6erqkpqbG/HTv3l2OP/74ii93dZkeLFe5yQwBBBBAAAG7Bc4880zRniwdm5WXl2cC\npT59+lRC+fXXX+WYY46JWa8D548++mizTsdr6QB6HRB/8MEHx2ynU0Ro4OVlogfLS33yRgAB\nBBBAIKACaWlppmY6Tio6aeBz0UUXyZQpU0yQ1aJFC+ncuXP0Juax9mxFJw2m9LTgcccdZ1br\n4HZNeooxOn3//femN2zgwIHRq11/TIDlOjkZIoAAAgggEHyBwsJCU8n77rtP3njjjZgK62nC\n4uJieeGFF+SSSy4xY6liNggv6FxaOhO8znOlr9fTiNp7dcMNN5hNr7rqKmnTpo08+uij8thj\nj8ns2bNlwoQJ8qc//ckEWLfffnvFXbq6TIDlKjeZIYAAAgggYIeADlo/8sgjZcyYMeaKwOha\nt2vXzsxfVVZWJvFOD+q2t956qxkMf/jhh0uvXr3ksMMOk3/84x+i47A06ZWCOt2DnnIcPHiw\n6D71KkU9Pfjyyy9Lo0aNzHZe/fL2BKVXtSZfBBBAAAEEEEiqgAY4esXf+vXrRa8WrJj0VOFJ\nJ50kBxxwQMWnzLKOv7rnnntk4cKFZsyWXmlYMWkeOpZLT0PqvFgFBQXSpEmTuD1iFV+b7GUC\nrGQLs38EEEAAAQQsFnBOFUYTfPbZZzJjxgxzS5zo9fEeR8+FFe95XadTOBx66KG7etqT9Zwi\n9ISdTBFAAAEEELBPQE8Xaq/V6aefbk4R/vGPfwwsAgFWYJuWiiGAAAIIIOAvAb0Z89atW82k\nonprm3hTKTRt2tSMpdpvv/38VfgES5NwgKWXQ9500027zEZH+usllyUlJbvchicQQAABBBBA\nwD4B7bnSewzqzZ2LioriAhxxxBHm1KHexLk2pz0ag7VmzRozgEwrqgPWdC6KZcuWVaq3DjJ7\n9913ZcmSJbJt2zapU6dOpW1YgQACCCCAAAIIBF1gjwIsvZO1zkURnZo1axa9GPNYL8uMN6gt\nZiMWEEAAAQQQQACBgArsUYA1aNAgM6+E3jPoo48+Ep2qvuJdr9VHz6VqYFXxho0BtaNaCCCA\nAAIIIIBAXIE9CrB0Mi+d8EuT3sl6zpw58pe//CXuDlmJAAIIIIAAAgjYLrBHAVY0UpAvqYyu\nJ48RQAABBBBAAIG9FUg4wNKMXnvtNXn44YfNqUK9WjAUClXKX2duJSFQUWDj8g6SWXFlEpd3\nZLyVxL2zawQQQACBeAIbNmyItzqp6+rVq5fU/Se684QDrH//+9+ivVh6haBeStm4cWNfTEmf\naMXZHgEEEEAAAQQQSJZAwgHWpEmTzD2Fvv32WznooIOSVS72iwACCCCAAAK1XODAvw5ypQY/\n9R3pSj6JZJLwRKMrVqww09sTXCXCzLYIIIAAAgggYJNAwgGWzqyqvVc61T0JAQQQQAABBBBA\noLJAwgGWzn/VpEkT+d///d/I7O6Vd8saBBBAAAEEEEDAXoGEx2DpRKONGjWSESNGyOOPPy46\no3tubm4lwZkzZ1ZaxwoEEEAAAQQQQMAGgYQDLJ1+Yfv27XLsscf63iclJUV0klQ/JS2TJr+V\ny09GNVkWPzn78f1Yk9bx9pWWlmZWp6amWveep+4iauCn/8F479GaXqfvdU16ZxOS3QIJvwP+\n/Oc/i/7UhqRv9Hi9a16W3Tno+q1cXpokM28/OWvb+6k8yXR39u18odAPG9vq7nzQZmZmivPY\ncQn6X+c4l5WVZV2g4QRWOpVReXl50Jua+u1GIOEAazf78t1TZWVlUlxc7KtyNWjQQPSA68Uk\nbL6CcKkwfnHWg+3OnTuta3f9sMnOzjbjNTdu3OhSq/sjGw0u9Gfbtm2+Ow4lWygnJ8cc5/Ri\nKJ2M2qaUl5dneu02b94sev9ePyU9DpHcEwh0gOUeIzkhgAACCCCAQG0R0N7Fzz//XD755BNp\n0aKFXHjhhebLYE2WP+EA65FHHpHHHnusyjIsXry4ym3YAAEEEEAAAQQQcFNA5/M85phjRHta\nO3XqJE888YTcfffdMmPGDKlfv36NFSXhAKthw4Zy8MEHxxRAT8UtWbLE3JuwsLBQevfuHfM8\nCwgggAACCCCAgB8EdAaEAw44QKZNm2aKs2XLFjMjgnYg3XPPPTVWxIQDrEsuuUT0J1765Zdf\n5PTTT5eioqJ4T7MOAQQQQAABBBDwVEDHyd12222RMugFODqJ+sKFCyPrauJBwhON7i7TVq1a\nmUJrBKi9WiQEEEAAAQQQQMBPArfeeqt07949UqRVq1bJxx9/LO3bt4+sq4kHNRpgaYH2228/\n0asnFixYUBPlYx8IIIAAAggggEBSBHRez169ekmbNm2kX79+NZpHjQZYeknu6NGjzeRyzZs3\nr9GCsjMEEEAAAQQQQKCmBH777Tfp2rWraA/We++9Z6YWqal9634SHoM1ZswYee655yqVQef7\n0CsH161bJ33D9yvU0fkkBBBAAAEEEEDAbwLLly+X0047TfLz8+XTTz8VnaOyplPCAdaOHTtE\nR9xXTDpzb7t27cwg94EDB1Z8mmUEEEAAAQQQQMBzgaVLl8rJJ58shx9+uLz88suSrAlYEw6w\n+vfvL/pDQgABBBBAAAEEapvA1VdfbS7E086gr7/+OlJ8nQOrbdu2keXqPkg4wHIy1Nt+6Kj7\n+fPnm9sBHHnkkaI/9erVczbhLwIIIIAAAggg4BsBnU7qnXfeMeU59dRTY8qlVxZOmTIlZl11\nFvYqwPrmm2/MOKtZs2ZVynv48OFyyy23VFrPCgQQQAABBBBAwEsBnU4qFAq5UoSEAyy9eW7P\nnj3NjWt11lOdN6Ju3bqyaNEief7550Xnl9Cbuw4aNMiVCpAJAggggAACCCDgN4GEAyy9ilCD\nrG+//Tbmljk6WKxHjx5y1VVXyVNPPUWA5beWpjwIIIAAAggg4JpAwvNgzZw5U0455ZSY4Cq6\ntH/+85/NJKN6CSQJAQQQQAABBBCwUSDhAEunY9CpGnaVnOe4Vc6uhFiPAAIIIIAAAkEXSDjA\n0hsifvLJJ/Lll19WstGBYw8++KA0bNjQ3DKn0gasQAABBBBAAAEELBBIeAzWFVdcITq4XU8T\nXnnllXLccceZmVB1kPtf//pXMzZLB7uTEEAAAQQQQAABWwUSDrB0xtPPP/9cLr/8cnn88cdj\n3AoLC2XUqFFy6aWXxqxnAQEEEEAAAQTsE/ip70j7Kv1/NU44wNLXNWnSxEzG9euvv8rcuXPN\n/QcPOOAAczdqnbKBhAACCCCAAAII2CywVwFWeXm5ueFzixYtpFu3bsbvtddekyFDhsjNN98s\nZ5xxhs2m1B0BBBBAAAEEwgJlm7q64pCW/y9X8kkkk4QHuZeWlsrRRx8tOh3DTz/9FMlLry78\n6quv5KyzzpIJEyZE1vMAAQQQQAABBBCwTSDhAEvvP/jDDz/I22+/Lddcc03E69xzzxW9Q3XX\nrl3lhhtuEO3lIiGAAAIIIIAAAjYKJBxgvfnmm3LyySebnqqKYHon6uuvv15WrVolCxcurPg0\nywgggAACCCCAgBUCCQdYqpKRkbFLHA2yNGVmZu5yG55AAAEEEEAAAQSCLJBwgNW5c2f56KOP\nzFQNFWH0tOCIESOkcePGTDRaEYdlBBBAAAEEELBGIOGrCE8//XRp3769mWj0D3/4gxx55JGS\nl5cny5Ytk8mTJ8u8efPkpZdesgaQiiKAAAIIIIAAAhUFEg6wdJ6rf/3rX+YqQh2PFX3F4H77\n7WeWe/XqVTEflhFAAAEEEEAAAWsEEg6wVCY7O1tefPFF0XsP6mB27b3af//9pWnTppKSkmIN\nHhVFAAEEEEAAgdonoFNOaWfRrFmz5IQTTpCOHTvWeCUSHoMVXQINplq1aiWdOnWSZs2aEVxF\n4/AYAQQQQAABBHwnsHbtWjNO/PbbbzcBls7fOXjw4Bov5171YNV4KdghAggggAACCCDggsDw\n4cPNWbfp06eb3KZMmSJnnnmmDBw4sEYv0KtWD5YLDmSBAAIIIIAAAgjUmMD5558vzz77bGR/\n++yzj3msc3jWZKIHqyY12RcCCCCAAAII+FrgxBNPNOUrKSmRTz/9VO6++24z1ElvA1iTyTc9\nWGVlZfLCCy/Ipk2bKtVv8+bN8t5778mkSZNkyZIllZ5nBQIIIIAAAgggkIjAmDFj5IILLpAZ\nM2bIoEGDJDW1ZkOimt1bIjWrsO3o0aNl7NixUlxcHPOMXqXYs2dPM8eWjva/7LLLDEbMRiwg\ngAACCCCAAAIJCFx33XWyfv16efnll02gpbMj1GTyPMDSc55DhgyRN954I2697rvvPunRo4do\npDls2DC5+OKLZeTIkWaKiLgvYCUCCCCAAAIIILAHAunp6XLhhRdKt27dzFmyPXjJHm/i+Ris\n+++/X9LS0uSBBx4wXXTRJV+3bp3MnTtXbrnllsgUEGeffbbp6ZozZ460bds2svnOnTtlw4YN\nkWV9oKcda7rLLyaDaiz4tVzVqJIvX+onZ53WxE/lcaPBnHnxqLvn32XdaO5IHrS7mP912/7f\nI28Anz847bTT5JxzzjFXDTpF3bhxozj3UnbWVfev5wHW0KFDRUfwL168uFJdVq5cadY1adIk\n8lyDBg3MjaRXr14dE2D9/PPPpqcrsmH4wZ133im9e/eOXuWbx85VC74pUEAL4idnvQG6n8rj\nZpPXqVNH9MfGpHe/0B8bU0FBgeiPjammP6xtNExWnfWsmHbunHrqqXLggQeaidN1ygYd512T\nyfMAa3cfOCtWrJCsrCzzE11pvfehnjeNTrpOo9LoVFRUJNu2bYte5flj/ZDVbzV+K5fnMEkq\ngF+c9e4H65Ycm6Raxt9tbuNp8Z9wca32ZOj/sPYw649NSf/P9f9dZ4zW3nSbkp6VyMjIkB07\ndkh5eblNVRc95aQ/27dv991QFj0OkUSuvPJK+eyzz+Twww83X/z0vfrkk0+acVg16eN5gLW7\nymil4x2U9WCVk5MT81Lt5Ro1alTMOj1lWDEQi9nAgwWnB85v5fKAwpUs/eKswb7byQ911w+a\nRo0amQ9a7YK3KWlgqb0YGuRXvHgn6A56fNaeq61bt4peCm9T0i/72mOpV79rcO2n5MVxyE/1\nd8qivemvvvqq6DFJhyK1aNHCDFVynq+pv74OsBo2bGi++ek/aXRApVM58EapqbcA+0EAAQQQ\nQMA+gWSfwvb1yEu9v6F+A549e3ak5XXQu3Y5R4/LijzJAwQQQAABBBBAwAcCvg6wNLrUSyfH\njRtnuti1q13nyurevbs57eADP4qAAAIIIIAAAghUEvB1gKWl7devnxkoqpdUnnvuuaZH69pr\nr61UEVYggAACCCCAAAJ+EfDNGCwdZDZtWuWrngoLC+XRRx81t9DRK1Nyc3P9Ykc5EEAAAQQQ\nQACBuAK+CbDili5qZX5+ftQSDxFAAAEEEEAAAf8K+P4UoX/pKBkCCCCAAAIIIBBfoNb0YMUv\nPmsRQAABBBBAwK8Cafn/8mvRkl4uAqykE5MBAggggAACdgnUq1fPrgrHqS2nCOOgsAoBBBBA\nAAEEEKiOAAFWdfR4LQIIIIAAAgggEEeAACsOCqsQQAABBBBAAIHqCBBgVUeP1yKAAAIIIIAA\nAnEECLDioLAKAQQQQAABBBCojgBXEVZHj9ciUAsFDpt4q+ulXvLnUa7nSYYIIICAlwL0YHmp\nT94IIIAAAgggEEgBAqxANiuVQgABBBBAAAEvBQiwvNQnbwQQQAABBBAIpAABViCblUohgAAC\nCCCAgJcCBFhe6pM3AggggAACCARSgAArkM1KpRBAAAEEEEDASwECLC/1yRsBBBBAAAEEAilA\ngBXIZqVSCCCAAAIIIOClAAGWl/rkjQACCCCAAAKBFCDACmSzUikEEEAAAQQQ8FKAAMtLffJG\nAAEEEEAAgUAKEGAFslmpFAIIIIAAAgh4KUCA5aU+eSOAAAIIIIBAIAUIsALZrFQKAQQQQAAB\nBLwUSPcyc/JGAAE7BErWniyZLlZ1R8ZbLuZGVggggEBlAXqwKpuwBgEEEEAAAQQQqJYAAVa1\n+HgxAggggAACCCBQWYAAq7IJaxBAAAEEEEAAgWoJEGBVi48XI4AAAggggAAClQUIsCqbsAYB\nBBBAAAEEEKiWAAFWtfh4MQIIIIAAAgggUFkg8NM0pKSkVK61D9b4tVw+oKnRItjsTN1r9K1U\nrZ3Z1hbR9Y1+XC3EWvJip77613lcS4pOMWtYINABVnp6ujRo0KCGyaq3u7S0NLMDv5WrerXy\n76ttdqbu3r4vnQ/XnJwcyc7O9rYwLueemvrfkyN169aV3Nxcl3P3Njun7gUFBRIKhbwtDLl7\nKhDoAGvnzp2yefNmT4ErZq4fepmZmbJ27dqKT7GcBAG/OBcVFSWhdrvfpV/qvvtSJudZP9Q9\nKytL6tevL1u3bpXi4uLkVNSne9WgUgMMrXdJSYlPS5mcYuXl5YkGlhs2bJDS0tLkZLKXe/Xi\nOLSXRQ3EyxiDFYhmpBIIIIAAAggg4CcBAiw/tQZlQQABBBBAAIFACBBgBaIZqQQCCCCAAAII\n+EmAAMtPrUFZEEAAAQQQQCAQAgRYgWhGKoEAAggggAACfhIgwPJTa1AWBBBAAAEEEAiEQKCn\naQhEC1GJpAkcNvHWpO073o7n/z7eWtYhgAACCARRgB6sILYqdUIAAQQQQAABTwUIsDzlJ3ME\nEEAAAQQQCKIAAVYQW5U6IYAAAggggICnAgRYnvKTOQIIIIAAAggEUYAAK4itSp0QQAABBBBA\nwFMBAixP+ckcAQQQQAABBIIoQIAVxFalTggggAACCCDgqQABlqf8ZI4AAggggAACQRQgwApi\nq1InBBBAAAEEEPBUgADLU34yRwABBBBAAIEgChBgBbFVqRMCCCCAAAIIeCpAgOUpP5kjgAAC\nCCCAQBAFCLCC2KrUCQEEEEAAAQQ8FSDA8pSfzBFAAAEEEEAgiAIEWEFsVeqEAAIIIIAAAp4K\nEGB5yk/mCCCAAAIIIBBEAQKsILYqdUIAAQQQQAABTwUIsDzlJ3MEEEAAAQQQCKIAAVYQW5U6\nIYAAAggggICnAgRYnvKTOQIIIIAAAggEUYAAK4itSp0QQAABBBBAwFMBAixP+ckcAQQQQAAB\nBIIoQIAVxFalTggggAACCCDgqQABlqf8ZI4AAggggAACQRQgwApiq1InBBBAAAEEEPBUgADL\nU34yRwABBBBAAIEgCqT7vVKbN2+W6dOnVypm586dJSMjo9J6ViCAAAIIIIAAAl4L+D7Amjlz\npgwfPlwaNmwYY9WhQwcCrBgRFhBAAAEEEEDALwK+D7AWLFggbdu2lVGjRvnFjHIggAACCCCA\nAAK7FfD9GCwNsFq3br3bSvAkAggggAACCCDgJ4Fa0YOVlZUlQ4cOlXnz5kmbNm1kwIAB0rRp\n0xjHn3/+WXr16hWzbvDgwXLBBRfErPN6ISUlxRRhn3328booVuRvszN198dbPDc3V/THxpSf\nny/6Y1NyjvH169e3qdrUNY6ArwMsHeC+cuVK2XfffU3w1LFjR5k8ebL0799fxo8fL3Xr1o1U\nKTU1VQoLCyPL+kADs/Ly8ph18RY2Lu8Qb3XS1jVo/tUelStpBbBox3vS/kHloO7etqzzQRsK\nhUR/bEo2110/i7T+fvz/07KR3BPwdYClAdSkSZNEvwlkZmYalUMPPVT69OkjU6dOlZ49e0ak\n9t9/f3n//fcjy/pgw4YNsmbNmph18Rb+u+d4zyRv3Z6UK3m527Nnm52pu7fvc/2Cp8eurVu3\nSnFxsbeFcTn3nJwcKSgoEP2SXFJS4nLu3maXl5dnvvzr509paam3hamQe1FRUYU1LCZTwNcB\nln4L0N6r6NSqVStp1KiRrFixIno1jxFAAAEEEEAAAd8I+Lq/cNGiRaa3aunSpREwDaz0m3nF\nMViRDXiAAAIIIIAAAgh4LODrAKtly5aSnZ0tTz/9tKxfv970Wo0ePdqMterSpYvHdGSPAAII\nIIAAAgjEF/B1gKVFHjRokCxcuFDOO+88M9B92bJl8uSTT4qe4ychgAACCCCAAAJ+FPD1GCwF\nO+SQQ2TChAmydu1aM3O7DpwkIYAAAggggAACfhbwfYDl4FW8VY6znr8IIIAAAggggIDfBHx/\nitBvYJQHAQQQQAABBBCoSoAAqyohnkcAAQQQQAABBBIUIMBKEIzNEUAAAQQQQACBqgQIsKoS\n4nkEEEAAAQQQQCBBAQKsBMHYHAEEEEAAAQQQqEqAAKsqIZ5HAAEEEEAAAQQSFCDAShCMzRFA\nAAEEEEAAgaoECLCqEuJ5BBBAAAEEEEAgQYFaM9FogvVicwQQQACBKIHM0h5RS8l9uHNjeP8F\nXyU3E/aOgM8F6MHyeQNRPAQQQAABBBCofQIEWLWvzSgxAggggAACCPhcgADL5w1E8RBAAAEE\nEECg9gkQYNW+NqPECCCAAAIIIOBzAQIsnzcQxUMAAQQQQACB2idAgFX72owSI4AAAggggIDP\nBQiwfN5AFA8BBBBAAAEEap8AAVbtazNKjAACCCCAAAI+FyDA8nkDUTwEEEAAAQQQqH0CBFi1\nr80oMQIIIIAAAgj4XIAAy+cNRPEQQAABBBBAoPYJEGDVvjajxAgggAACCCDgcwECLJ83EMVD\nAAEEEEAAgdonQIBV+9qMEiOAAAIIIICAzwUIsHzeQBQPAQQQQAABBGqfAAFW7WszSowAAggg\ngAACPhcgwPJ5A1E8BBBAAAEEEKh9AgRYta/NKDECCCCAAAII+FyAAMvnDUTxEEAAAQQQQKD2\nCRBg1b42o8QIIIAAAggg4HMBAiyfNxDFQwABBBBAAIHaJ0CAVfvajBIjgAACCCCAgM8F0n1e\nvmoVLyMjQ/Ly8qrcx4ZlVW5S4xs0bty4xvfJDisL2OxM3Su/H7xYk5ubKzk5OV5kHZOnF8e5\n/Pz8PToGxxS0li+kpKSYGhQWFtbymlD86goEOsAqLS2VTZs2VWmUWeUWNb/B6tWra36n7LGS\ngM3O1L3S28HVFVlZWVK/fn3ZsmWLFBcXu5p3vMy8OM7p8bekpCRecQK7Tr/U161bV9avXy/6\nGeSnVFRU5KfiBL4snCIMfBNTQQQQQAABBBBwW4AAy21x8kMAAQQQQACBwAsQYAW+iakgAggg\ngAACCLgtQIDltjj5IYAAAggggEDgBQiwAt/EVBABBBBAAAEE3BYgwHJbnPwQQAABBBBAIPAC\nBFiBb2IqiAACCCCAAAJuCxBguS1OfggggAACCCAQeAECrMA3MRVEAAEEEEAAAbcFCLDcFic/\nBBBAAAEEEAi8AAFW4JuYCiKAAAIIIICA2wIEWG6Lkx8CCCCAAAIIBF6AACvwTUwFEUAAAQQQ\nQMBtAQIst8XJDwEEEEAAAQQCL0CAFfgmpoIIIIAAAggg4LYAAZbb4uSHAAIIIIAAAoEXIMAK\nfBNTQQQQQAABBBBwW4AAy21x8kMAAQQQQACBwAsQYAW+iakgAggggAACCLgtQIDltjj5IYAA\nAggggEDgBQiwAt/EVBABBBBAAAEE3BYgwHJbnPwQQAABBBBAIPACBFiBb2IqiAACCCCAAAJu\nCxBguS1OfggggAACCCAQeAECrMA3MRVEAAEEEEAAAbcFCLDcFic/BBBAAAEEEAi8AAFW4JuY\nCiKAAAIIIICA2wIEWG6Lkx8CCCCAAAIIBF6AACvwTUwFEUAAAQQQQMBtAQIst8XJDwEEEEAA\nAQQCL0CAFfgmpoIIIIAAAggg4LYAAZbb4uSHAAIIIIAAAoEXIMAKfBNTQQQQQAABBBBwWyDd\n7Qz3Jr/NmzfL559/Lvq3ffv20rx5873ZDa9BAAEEEEAAAQRcEfB9D9bChQulZ8+eMnnyZJk1\na5ZcdtllMmPGDFdwyAQBBBBAAAEEENgbAd/3YN13333So0cPGThwoKSkpMgLL7wgI0eOlFde\necUs702leQ0CCCCAAAIIIJBMAV/3YK1bt07mzp1rerA0uNJ09tlny/Lly2XOnDnJdGHfCCCA\nAAIIIIDAXgv4ugdr5cqVpmJNmjSJVLBBgwaSmZkpq1evlrZt20bWL126VG6++ebIsj64+OKL\npUuXLjHr4i1s/m828Z5K2jqtByn5AjY7U/fkv792l4PzpbBOnTqSlZW1u01dec6L41zdunUl\nJyfHlfr5JZO0tDRTlIKCAgmFQn4pFuXwQMDXAdaKFSvMganiwSkvL0/Wr18fw7V161b55ptv\nYtadddZZJhiLWRlnoUHzr+KsTe4qDRL9kNZeP84PxfCkDLbW3dZ6e/ImC2eanu6Pw6wXxzm/\n1N2Lts/IyPAiW/L0kYA//vN3AaJv0J07d1Z6tqysrNK3otatW8vs2bNjttWrDjVI81NyeuD8\nVi43jBo3bmy+0a1Zs8aN7HyVR1FRkWzfvl1+++03X5Ur2YXRD9hGjRqJfgHauHFjsrPz1f71\ni2H9+vXN1c/FxcW+KluyC6O9VtqDs2HDBikpKUl2dr7av3YAaM/d2rVrpbS01Fdl0+MQyT0B\nXwdYDRs2FA2m9OAc3c28adMmifdGqfhtyemid4+TnBBAAAEEEEAAARFfD3Jv1qyZ6V6P7pnS\nQe/l5eUSPS6LhkQAAQQQQAABBPwk4OsAS7uYu3XrJuPGjRPtYt+2bZuMHTtWunfvbk47+AmS\nsiCAAAIIIIAAAo6ArwMsLWS/fv3MQPVzzjlHzj33XNOjde211zrl5y8CCCCAAAIIIOA7AV+P\nwVKtwsJCefTRR0XHXenlr7m5ub5DpEAIIIAAAggggEC0gO8DLKew+fn5zkP+IoAAAggggAAC\nvhbw/SlCX+tROAQQQAABBBBAII4AAVYcFFYhgAACCCCAAALVESDAqo4er0UAAQQQQAABBOII\nEGDFQWEVAggggAACCCBQHQECrOro8VoEEEAAAQQQQCCOAAFWHBRWIYAAAggggAAC1REgwKqO\nHq9FAAEEEEAAAQTiCBBgxUFhFQIIIIAAAgggUB0BAqzq6PFaBBBAAAEEEEAgjgABVhwUViGA\nAAIIIIAAAtURIMCqjh6vRQABBBBAAAEE4gjUmnsRxil7lavq1asn+uOn1L9/f5k7d658+OGH\nfiqWK2W54IILJDMzUyZMmOBKfn7JZMeOHXLqqafKkUceKY888ohfiuVKOZYsWWLqfsYZZ8iQ\nIUNcydMvmXzxxRdyyy23yGWXXSYXXXSRX4rlSjneeOMNefzxx039u3bt6kqefsnkmWeekYkT\nJ5r6t2vXzi/FohweCAQ6wPLAs8os16xZI8uWLatyuyBusGLFCsnKygpi1XZbp1AoZNq8WbNm\nu90uiE+Wlpaauq9fvz6I1dttnbZt22bqvmnTpt1uF8Qni4uLTd23bt0axOrttk4bNmwwdd++\nfftut+PJ4AtwijD4bUwNEUAAAQQQQMBlAQIsl8HJDgEEEEAAAQSCL8ApQpfb+He/+500aNDA\n5Vz9kV3Hjh0lIyPDH4VxsRSpqalmHFLr1q1dzNUfWeXm5pq6H3roof4okIul0P9zHXvXsmVL\nF3P1R1b77befqXtRUZE/CuRiKQ466CBTd7+N/3WRgKz+TyAlPD4khAYCCCCAAAIIIIBAzQlw\nirDmLNkTAggggAACCCBgBAiweCMggAACCCCAAAI1LMAYrBoG3dXutmzZIlOnTpWffvpJ0tLS\n5JhjjpGjjjpKdIxKkNNLL70k0Zcr67iE/fffX4444gjRsUk2JJ33bPr06XL66adL06ZNbaiy\nvPXWW7J27dq4de3WrZsEfcqKhQsXyg8//CBz5swxc/HpuJxTTjnF/O/HRQnAyn//+9+mzXv0\n6BFTm1WrVsk777wjBx54oJx00kkxzwVh4eWXX5aSkhLp3r27NGnSpFKVXnvtNdm4caP86U9/\nkpycnErPsyK4AgRYLrTt999/L3fddZfocDcd5F5WViYPPvigCTx0Mr4gD37WSUXr1q0r++67\nr5SXl5sDza+//irt27c3JjYMen/++edl5syZovPj3HDDDS6847zPQgOsdevWSfPmzSsV5vjj\nj6+0LkgrHnvsMXn99dfl6KOPlhYtWsjy5ctNwDl58mS59957pX79+kGqbqQuGmDNmzdPogMs\nDa6uvfZaKSgokAsvvDCybZAe6DFO/7f1uH7llVfGVG316tXy6KOPmnXnnHMOAVaMTvAXCLCS\n3MabN2+W2267TY499li5+eabIxNtarClQdfgwYNl9OjRolfdBDXpAbd3796R6ulBWA9EX3/9\ntXTo0CGyPogPtBfnyy+/FJ3B/7nnnpN+/fpZc5A98cQTzfs7iO26qzppMP3BBx+IzuYd/cVJ\nP2g10Pjb3/4mAwcO3NXLA7XeCa4aNmwoI0aMCHRvvR6/P/roo0oBlt6xQ+u/q97cQDU4lakk\nYMc5mkrVdm+FniJLSUmRG2+8MRJcae66Tm+jkZeXZ26r4F6JvM9JP3jq1KljxUFnypQp5rTg\neeedZ04P6YcvKZgCOmv5Cy+8YILo6OBKa9u4cWOz3pY7GaxcuVIGDBhgeq4ffvjhQAdX2r56\n+ld75nUISHTS//cuXbpEr+KxRQIEWElu7NmzZ+9yrJXel09PI+g4jSAnHZ+gt0r57bffZPHi\nxTJmzBhzwD3hhBOCXG1TNx17ovdi01OhOieSnjqyJWkvrZ4Wjv4J8qwwOtZO66q91fFS586d\nTZAV77kgrdPgSnvrnJ4r/TIV9KSnffVeo9qL5SS9JZoGXUE/Je7Ul7+VBThFWNmkRtcsXbpU\n2rZtu8t96sDPt99+W/SGwBpwBTHpt3r9iU56ujToE67+5z//Mfck08Gvms4880x58803RYPu\n3b0nop1q82Mdh6U/0emSSy6pdBol+vna/Fi/KOl4Q+2tcpKeGvruu++cRfNXv1gE9eIWHYuk\nwZXeh1G/TOng7miPGIiALWhPld7k2RmHpb1XnTp1CuxxPWDNl5TqEGAlhfX/71Svmou+iu7/\nP/PfR9q7o6cJgxpcaS316hkd4KlJx6TpB5EOBNYbwupzQU3ae6VXFelpA+fUgbb1G2+8YUWA\npR8uvXr1imneIH/Y6v+6Xi2spwqdq8WWLFkiTz/9tDHYuXOn6cUdP358YAMsHXelPXVDhgyR\nK664QoYNGyZ6IY9eOR30dPLJJ8vIkSPN/7p+cdarxnXsJcleAQKsJLf9YYcdFvlwjZeVfvC2\na9cu3lOBWacfPNFXk2nvjXafa29OUAMs/ZDV0wXaSzd27NhIW2rPhQ581W/5+fn5kfVBfKCn\nTfT9b0vSuuop0Pnz55thAVpvHQKgl+lr0vUadAQ56WDvv/zlLyaguuOOO0wMKU+7AAAMtElE\nQVSAMW7cuMDXW9tUj3Pa3vp/rwGlDovQ6XiCPgQkyO/n6taNMVjVFazi9Xr/PT1VpFfMVUw6\nV45+yznuuOMqPmXFcmlpaWDrqe2qH7Z65aBzilT/6tVlejn3u+++G9i621oxnZJBeyy1zbVn\numIK8vvdqav23Dm9VfrFUU8Jv/jii/Ltt986mwT6r54m1ABLv0RpT156On0YgW7wKipHgFUF\nUHWf1vmerr76anPF4N///nfRLnSdH+if//ynXHfddXLWWWeJXmEW5LRixQozD5TOBfXNN9+I\nniLR3isd/B3UpOPqdJoCHZMTnbRXR6em0LFJQR7wHV1nWx7rlcF6ikjf74MGDZL333/fDHKe\nNWuWec/fdNNNZjySnia2JfXp00fatGljpqTRHp2gJ51IVec904tZuHow6K1ddf0Ir6s2qvYW\n//M//yPZ2dnmW42OR9AejAMOOMAEVpdeemm19+/3HWgwpT+a1EHH4eg324svvtjvRd+r8i1a\ntMicFtA6xksaVOsUHRps6sSzpOAIaA/WE088IToflvZkabClqaioSC644AIzH5wtUzVovbUH\nR08V6nFOJ1nV+bA0EA1q0uBZz0j8/PPPcvjhhwe1mtRrDwVSwt+iQ3u4LZvVgICOzdHBrkEf\nf1MDVOwCgVovoFcR6jQFQb1qsNY3EBVAIIkCBFhJxGXXCCCAAAIIIGCnAGOw7Gx3ao0AAggg\ngAACSRQgwEoiLrtGAAEEEEAAATsFCLDsbHdqjQACCCCAAAJJFCDASiIuu0YAAQQQQAABOwUI\nsOxsd2qNAAIIIIAAAkkUIMBKIi67RgCBWAGdpkRvAqw3A05G0kkeV65cmYxds08EEEAgIQEC\nrIS42BgBewT01i4PPvigrFmzpsYqPWXKFGnZsqV8/PHHNbbP6B2dfvrpcu6550av4jECCCDg\niQABlifsZIqA/wV01u2bb75ZduzY4f/CUkIEEEDAZwIEWD5rEIqDgF8E9I4DJAQQQACBvRPg\nXoR758arEKgxgaVLl8rYsWPN/QsLCgrksMMOk8svvzxyo+gFCxaYmwVfc8015n6WeqPwwsJC\nc6PwU045RdLS0mLKooHRuHHj5MsvvxQd83TUUUfJlVdeKbrv6FRcXGxuSPzRRx+Z1T179hTd\nX0ZGhkycONHkpU889NBDZh96b8WnnnrK3Etyv/32k1GjRpkb+WpZGzVqZPbx7rvvyrRp00TL\nXK9ePWnbtq3Ju+JNr83Ge/Frd2Xe1e70FOeLL74o8+bNE73hsN4H9Oyzz5ZOnTpFXlJeXm5u\n0Pvee+/Jb7/9JgcffLCceeaZMdvoxlW1VWSHPEAAAQT0XoQkBBDwRiAciIQaNGgQCgcjofAH\neij8oR8K3yA3FA4CQuH72JlCvfPOO3q/0NAZZ5xhtu3bt2+oa9euZt2gQYNiCr569epQ+AbS\n5rlwkBAKj0cy+27RokVo9uzZkW1LSkpCHTt2DIWDqVC3bt1CnTt3DqWmpoaOOOKIUDhAC915\n552h8Fgps58TTjghNGTIEPPaI488MhQOwkLhgMpsr6/59ddfzXPhm5pH8v39738f2nfffc3y\nQQcdFNq+fbvZZvLkyWZdeCxWpCx7+qCqMut+2rVrF2rfvn1kl5999pkpazjAM2Za53BAGgrf\ncDgUvhlzZLvrrrvOrAsHt6HwTZlD4Zs2m+VHHnkkss2etFVkYx4ggID1AmK9AAAIeCiggUtm\nZmYo3DMSKcWrr75qgpAnnnjCrHMCLA1qfvnll8h2w4cPN9v9/e9/j6y77LLLKq0LX7UXKioq\nMsGbs+HFF19sAjkNQJz09ttvm9c+88wzZtWwYcPMshNA6UoNsDTYu+GGG0Lh3rFQuFfIbPvh\nhx+a9TfddJNZ1l/hXqHQ1Vdfbda/9dZbZn11Aqw9KXPFAOukk04K5eXlhcJXFkbKpdYaxB59\n9NFm3aZNm0zQFe4hjGxTVlYW0mBrn332MQGnPrEnbRXZAQ8QQMB6AcZg0YmJgIcC4R4gM4g8\nHOhI+EPdlOTCCy+UFStWyIABA2JKdv3118v+++8fWXfjjTdKw4YN5bXXXjPrNmzYYE4NdujQ\nQc4777zIds2bN5dw75I5dff999/rlyp588035Q9/+IOceOKJke3OOussc9pPr/LbXdIy33vv\nvVKnTh1p3bq12VTLNWHCBLntttsiLw33Ekm4J8ssV/dKxL0ps77m1ltvlX/84x8SDpQi5WrW\nrJmEe7liro7UOoV7+Iy7bqjLn3zyiSxatChyCjaRtopkxgMEELBWgDFY1jY9FfeDgI5fGj9+\nvPTq1Uv69+8v4VN/ZmyVjoeqmDQoiE7hni855JBD5LvvvjOrddyTBhXhHhkTPEVvG+6FMos/\n/vijGdul24R7o6I3MY91nFdVScdfZWdnx2ymQZn+fPXVV/LNN9/I3Llzzc+MGTPMdtW9EnHh\nwoWmXomUWQM8nbZBx11pEDpnzhyZP3++8dLHTtAV7uGS8KlWMyVF06ZNJdyzJeHTsWa6h2OO\nOSZSz0TaKvIiHiCAgLUC9GBZ2/RU3A8C4fFJJiC55557RB9PmjRJdDC5Pv78889jiqiDxium\nnJwc2bhxo1kdHrNl/mrPkva2RP9oL9Yf//hH0WBi2bJlZjt9vDcpPGas0ss0YAufjpPjjjtO\nwqcPTaB14IEHyuDBgyttuzcr9rbMf/vb30QDwvC4Knn++efNoP8+ffpIdOCk5XnggQdEB+hr\nT9+SJUtE2yM8lk10W+dqykTaam/qyGsQQCBgAuFvvCQEEPBQwBkArkVYt25d6OmnnzYDyHVA\ntiZnDNYbb7xhlqN/hU/Rhc455xyzSsdDhQ9PIR1sXjHpwHUn6Rgk3e722293VkX+vv766yEd\nJ6VpV2OwnLFLkReFH4RPZ5p9jhkzJhTurYo8FQ4YzfrRo0ebdXs7BmtPyxw9BkvHXelYKx30\nr+PQolP46sZQ48aNI6vUxzHSsWPhKzBDOn5LnT744IPIdlW1VWRDHiCAgPUC9GAFLGCmOrVL\nQHuVtIdly5YtpuD169eXq666ypz6q3jLl1deeSWmcjoNg57y0l4jTa1atZLwlXtmugHtUYpO\nvXv3NtMm6G1q9DSY5hkOdkSnJ3CSTk+g24UH15tVzvQPe3J6T08Nam+a9vjoNA9OCgeH5qHT\nC+SsT/TvnpY5er966lTz1fFo2oPnJD2VqqdKnTL95z//kdzcXBk5cqTZRE8tHnvssZExcE47\nJNJWTl78RQABewUIsOxte2ruAwEdqB6eWkEuuugi0WDkiy++MLOn6xghHewenXRuKp1Z/Ycf\nfpBwb5acf/75ZryQnpLTpIGNzr4ens7AjB/SQdoa+Gge+trwVAQSnq5BNIDQU2I6L5QOdNc8\n9RY2Og5MX6t5aNK5tjTdd999Jj+zsItfOjZK59y65ZZbzBxYOvZKx5S9/PLL5hXOacxdvLzK\n1Xta5ugd6RxcahKejsHUT4NLDVJPO+00c/p08+bNZsyalr1Lly5y//33myBr5syZptx33XWX\nGa+m47E0JdJW0eXgMQIIWCpgfR8eAAh4LBAe7xPSeZrChyDzk5+fHwpf/RbSqQI0OacIw1fo\nhcK9VGab8CBzMy9WeGxSpdKHgykzj5OzPz1NptM3RJ/e0he99NJL5jSZs51OSaDrnKRzaoWD\nD5OfzqOlSZfjnSLUU5tXXHFFZH8611R4Ms9QeHC6mSLi1FNPNa/f21OE5sXhX1WVOfoUob5G\np7wIXylpTrlqPfW04LPPPhsKT5hq6vXpp5+aXYd7As38WY5FOKAzdZ01a5Z53vlVVVs52/EX\nAQQQSFGC8EGFhAACHgroabjwOCNTAj3Vpz02TtLB1zqFQni+K3O6S6+oCwcK5rSWs028v3pq\nKxz4SMvw1X16CmxXSfPV/HWqBR0YXzHpVXh61aAOnq8q6SlHPf2mPWV7sn1V+9vV81WVueLr\nVq1aJXrzap2iYXdJLxTQKTLUbFcXAeyurXa3b55DAAG7BAiw7GpvalsLBSoGWLWwChQZAQQQ\nsE6AebCsa3IqjIA/BMIzxpspEaoqjfYkDR06tKrNeB4BBBDwlQABlq+ag8IgUFlAr6DTq/v0\nyr8gJT2FqQPPq0rx5v+q6jU8jwACCHgtwClCr1uA/BFAAAEEEEAgcAKVR7QGropUCAEEEEAA\nAQQQcFeAAMtdb3JDAAEEEEAAAQsECLAsaGSqiAACCCCAAALuChBguetNbggggAACCCBggQAB\nlgWNTBURQAABBBBAwF0BAix3vckNAQQQQAABBCwQIMCyoJGpIgIIIIAAAgi4K0CA5a43uSGA\nAAIIIICABQL/D5JQtNHN2qoVAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "data %>% \n",
    "    filter(type %in% c(2,3)) %>% \n",
    "    select(type,spectral_class) %>%\n",
    "    ggplot(aes(x=spectral_class, fill=type)) +\n",
    "        geom_bar(position = position_dodge(preserve = \"single\")) +\n",
    "        scale_fill_manual(values=cbPal[c(3,4)]) +\n",
    "        scale_x_discrete(drop=FALSE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: do types 2 and 3 share the same distribution of spectral class?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test differences in distributions over a *categorical variable* like spectral class, we can use [*Pearson's chi-squared test*](https://en.wikipedia.org/wiki/Chi-squared_test):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chi-squared test for statistical independence\n",
    "\n",
    "#### Theory\n",
    "\n",
    "$H_0$: Probability distribution for the categorical variable (e.g. spectral class) is independent of group (e.g. star type).\n",
    "\n",
    "$H_1$: Probability distribution depends on group.\n",
    "\n",
    "We need to find out whether the differences in observed frequencies between the two groups are small enough to have arisen by chance. We do this by constructing a [*contingency table*](https://en.wikipedia.org/wiki/Contingency_table) showing the observed frequency of each outcome for each of the two groups, and comparing to the *expected frequencies* under $H_0$.\n",
    "\n",
    "\n",
    "The test statistic is\n",
    "\n",
    "$$X^2 = \\sum^k_{i=1}{\\frac{(x_i-m_i)^2}{m_i}}$$,\n",
    "\n",
    "where $k$ is the number of cells in the table and $x_i$ and $m_i$ are the observed and expected frequencies for each cell.\n",
    "\n",
    "Under $H_0$, $X^2$ follows a [$\\chi^2$ distribution](https://en.wikipedia.org/wiki/Chi-squared_distribution), which is parametrised by the number of degrees of freedom. A contingency table of size $a \\times b$ has $(a-1)(b-1)$ degrees of freedom, i.e. the number of independent counts when row and column sums are held fixed.\n",
    "\n",
    "#### Assumptions\n",
    "\n",
    "- Sampling is random\n",
    "- Each sample is independent of the others\n",
    "- Expected frequency for each cell must be sufficiently large\n",
    "\n",
    "The approximation to the $\\chi^2$ distribution breaks down if expected frequencies are too low. It will normally be acceptable so long as no more than 20% of the events have expected frequencies below 5. Where there is only 1 degree of freedom (i.e. a 2 $\\times$ 2 table), the approximation is not reliable if expected frequencies fall below 10. \n",
    "\n",
    "\n",
    "#### Application"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will set $\\alpha=0.05$.\n",
    "\n",
    "We begin by making a contingency table for the observed spectral classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"OBSERVED:\"\n",
      "    spectral_class\n",
      "type  O  B  A  F  K\n",
      "   2  0 24  7  9  0\n",
      "   3  5 13 10  8  4\n"
     ]
    }
   ],
   "source": [
    "# OBSERVED\n",
    "\n",
    "obs <-\n",
    "    data %>% \n",
    "    filter(type %in% c(2,3)) %>% \n",
    "    select(type,spectral_class) %>%\n",
    "    table(exclude=list(\"G\",\"M\",0,1,4,5))\n",
    "\n",
    "print(\"OBSERVED:\")\n",
    "\n",
    "print(obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"row totals:\"\n",
      "type\n",
      " 2  3 \n",
      "40 40 \n",
      "[1] \"column totals:\"\n",
      "spectral_class\n",
      " O  B  A  F  K \n",
      " 5 37 17 17  4 \n",
      "[1] \"total observations\"\n",
      "[1] 80\n"
     ]
    }
   ],
   "source": [
    "print(\"row totals:\")\n",
    "tot_rows <- margin.table(obs,1)\n",
    "show(tot_rows)\n",
    "\n",
    "print(\"column totals:\")\n",
    "tot_cols <- margin.table(obs,2)\n",
    "print(tot_cols)\n",
    "\n",
    "print(\"total observations\")\n",
    "tot_obs <- sum(tot_rows)\n",
    "print(tot_obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the row and column totals, we can calculate the *expected* frequencies under $H_0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"EXPECTED:\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>O</th><th scope=col>B</th><th scope=col>A</th><th scope=col>F</th><th scope=col>K</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>2</th><td>2.5 </td><td>18.5</td><td>8.5 </td><td>8.5 </td><td>2   </td></tr>\n",
       "\t<tr><th scope=row>3</th><td>2.5 </td><td>18.5</td><td>8.5 </td><td>8.5 </td><td>2   </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllll}\n",
       "  & O & B & A & F & K\\\\\n",
       "\\hline\n",
       "\t2 & 2.5  & 18.5 & 8.5  & 8.5  & 2   \\\\\n",
       "\t3 & 2.5  & 18.5 & 8.5  & 8.5  & 2   \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | O | B | A | F | K |\n",
       "|---|---|---|---|---|---|\n",
       "| 2 | 2.5  | 18.5 | 8.5  | 8.5  | 2    |\n",
       "| 3 | 2.5  | 18.5 | 8.5  | 8.5  | 2    |\n",
       "\n"
      ],
      "text/plain": [
       "  O   B    A   F   K\n",
       "2 2.5 18.5 8.5 8.5 2\n",
       "3 2.5 18.5 8.5 8.5 2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"EXPECTED:\")\n",
    "exp <- as.matrix(tot_rows) %*% t(as.matrix(tot_cols)) / tot_obs\n",
    "exp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected values for classes K and O are too small. We can combine columns O, F and K to get all of the expected values over 5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"OBSERVED:\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>A</th><th scope=col>B</th><th scope=col>OFK</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>2</th><td> 7</td><td>24</td><td> 9</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>10</td><td>13</td><td>17</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       "  & A & B & OFK\\\\\n",
       "\\hline\n",
       "\t2 &  7 & 24 &  9\\\\\n",
       "\t3 & 10 & 13 & 17\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | A | B | OFK |\n",
       "|---|---|---|---|\n",
       "| 2 |  7 | 24 |  9 |\n",
       "| 3 | 10 | 13 | 17 |\n",
       "\n"
      ],
      "text/plain": [
       "  A  B  OFK\n",
       "2  7 24  9 \n",
       "3 10 13 17 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# OBSERVED (combining columns)\n",
    "\n",
    "AB <- obs[,c('A','B')]\n",
    "OFK <- margin.table(obs[,c('O','F','K')],1)\n",
    "new_obs <- cbind(AB,OFK)\n",
    "print(\"OBSERVED:\")\n",
    "new_obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"EXPECTED:\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>A</th><th scope=col>B</th><th scope=col>OFK</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>2</th><td>8.5 </td><td>18.5</td><td>13  </td></tr>\n",
       "\t<tr><th scope=row>3</th><td>8.5 </td><td>18.5</td><td>13  </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       "  & A & B & OFK\\\\\n",
       "\\hline\n",
       "\t2 & 8.5  & 18.5 & 13  \\\\\n",
       "\t3 & 8.5  & 18.5 & 13  \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | A | B | OFK |\n",
       "|---|---|---|---|\n",
       "| 2 | 8.5  | 18.5 | 13   |\n",
       "| 3 | 8.5  | 18.5 | 13   |\n",
       "\n"
      ],
      "text/plain": [
       "  A   B    OFK\n",
       "2 8.5 18.5 13 \n",
       "3 8.5 18.5 13 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# EXPECTED (combining columns)\n",
    "\n",
    "AB <- exp[,c('A','B')]\n",
    "OFK <- margin.table(exp[,c('O','F','K')],1)\n",
    "new_exp <- cbind(AB,OFK)\n",
    "print(\"EXPECTED:\")\n",
    "new_exp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In R, the most convenient form of the chi-squared test for a contingency table is `chi2.test`, which just takes the table of observations as input and calculates the expected values and degrees of freedom accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tPearson's Chi-squared test\n",
       "\n",
       "data:  new_obs\n",
       "X-squared = 6.2612, df = 2, p-value = 0.04369\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "chisq.test(new_obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the chi-squared test, $p<\\alpha$, so there is enough evidence to reject $H_0$:\n",
    "it appears that white dwarves and main sequence stars follow different distributions for spectral class.\n",
    "\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
